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ABSTRACT 

Context. In addition to gas phase reactions, the chemical processes on the surfaces of interstellar dust grains are important for 
the energy and material budget of the interstellar medium. The stochasticity of these processes requires special care in modeling. 
Previously methods based on the modified rate equation, the master equation, the moment equation, and Monte Carlo simulations 
have been used. 

Aims. We attempt to develop a systematic and efficient way to model the gas-grain chemistry with a large reaction network as 
accurately as possible. 

Mettiods. We present a hybrid moment equation approach which is a general and automatic method where the generating function is 
used to generate the moment equations. For large reaction networks, the moment equation is cut off at the second order, and a switch 
scheme is used when the average population of certain species reaches 1. For small networks, the third order moments can also be 
utilized to achieve a higher accuracy. 

Results. For physical conditions in which the surface reactions are important, our method provides a major improvement over the 
rate equation approach, when benchmarked against the rigorous Monte Carlo results. For either very low or very high temperatures, 
or large grain radii, results from the rate equation are similar to those from our new approach. Our method is faster than the Monte 
Carlo approach, but slower than the rate equation approach. 

Conclusions. The hybrid moment equation approach with a cutoff and switch scheme is a very powerful way to solve gas-grain 
chemistry. It is applicable to large gas-grain networks, and is demonstrated to have a degree of accuracy high enough to be used for 
astrochemistry studies. Further work should be done to investigate how to cut off the hybrid moment equation selectively to make the 
computation faster, more accurate, and more stable, how to make the switch to rate equation more mathematically sound, and how to 
make the errors controllable. The layered structure of the grain mantle could also be incorporated into this approach, although a full 
implementation of the grain micro-physics appears to be difficult. 
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1. Introduction 

The chemistry of the interstellar medium can be roughly divided 
into two types: gas phase chemistry and grain surface chem- 
istry. The two types of chemistry are coupled by the adsorp- 
tion and desorption processes. Species adsorbed on the grain 
surface migrate in a random walk manner, and they may react 

. with each other upon encounter at the same site (a local potential 
minimum). The products can be released back to the gas phase 
through certain desorption mechanisms. In addition to the gas 
phase chemistry, grain chemistry is important for the material 

' and energy budget of the interstellar medium. For example, be- 
sides H2, molecules s uch as methanol are believed to be formed 
on the grain surfaces (iGarro d et al. 2007). because its relatively 
high abundance (see, e.g.. lMenten et al..,198&) cannot be repro- 
duced by gas phase chemistry. 

Several methods have been used to model the gas-grain 
chemistry. In the ra te equation (RE) approach (see, e.g., 
iHasegawa et alJl992h . the surface processes are treated the same 
way as the gas phase processes. This works fine when the num- 
ber of reactants on a single grain is large (under the assumption 



that the system is well stirred; see iGille spiel (l2007b ). but might 
not be accurate enough when the average population^ of some 
reactants on a single grain is small. This failure of the rate equa- 
tion is related to the treatment of two-body reaction. For the REs 
to be applicable, the probability of one reactant being present 
should be independent of that of another being present. This is 
not always true, especially when the average populations of both 
reactants are low, in which case they might be highly correlated. 
The flaws in em ploying the RE f or gra in-sur face chemistry were 
pointed out by ICharnlev et al.l (1 19971) and iTielens & ChamlevI 
(1997). 

To remedy this problem, modification schemes based on 
some empirical, heuristic, a nd/or physical reasoning have been 
applied to the RE approach (ICaselli et al.ll9 98: Stan tcheva et al.l 
2001), and are called modified rate equation (MRE) approach. 
The validity of this method has been questioned ( Rae et aT] 
2003). A modification scheme developed b v iGarrodi ( 2008 ) uses 
different functional forms for different surface populations, tak- 
ing various competition processes and refinements into account. 



* Member of the International Max Planck Research School (IMPRS) 
for Astronomy and Astrophysics at the Universities of Bonn and 
Cologne. 



Here by "population" we mean the number of a species in a volume 
of interest, and by "average" we mean an ensemble average (i.e. aver- 
age over many different realizations of the same system setup). Hence 
"population" can only take non-negative integer values, while "average 
population" is a non-negative real number. 
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It has bee n shown to work ver y well, even for very large reaction 
networks (iGarrod et al.ll2009l) . 

Mathematically, the gas-grain system should be vi ewed as a 
stochastic chemical sy stem (see, e.g., McOuanie 1967t lGillespiel 
1 1 9761: Iciiarnleyl 19981) . being described by a probability distribu- 
tion P(x, t), which is the probability that the system has a pop- 
ulation vector X at time t, with xi being the number of the /th 
species in the system. The evolution equation of P{x, t) is the 
so-called master equation, whose form is determined by the re- 
action network. 

Many sophisticated methods have been proposed (mainly 
outside the astro-ch emical community; see, e.g., the operator 



method described in Mattis & Glasserl (1 19981) . or the variational 



approach used bv lOhkubol (120081) 7 to solve the master equation 



However, these methods work fine only when either the chemi- 
cal network is small or some special assumptions are made in the 
derivation, thus their validity in the general case should be ques- 
tioned. It is unclear whether these methods can be generalized to 
large complex networks. 

The numerical so lution of the master equation has also 
been perform ed (Bihamet all |200lt (Stantcheva at al.J (20021 
IStantcheva & Herbst 2004). To limit the number of variables in 
the set of differential equations and to separate the deterministic 
and stochastic species, usually a priori knowledge of the sys- 
tem is required in these studies. The steady state solution of the 
master equation can also be obtained analytically in some very 
simple c ases, such as the formation of H? molecules on the grain 
surface (iGreen et aLlbOOUlBiham & Lipshtatll2002h . 

On the other hand, the master equation prescription can 
be "realized" through a stochastic simula tion algor it hm (S SA), 
proposed by Gillespid (11976.) (see also IGillespid ( 119771) and 
iGillespid (12007 )). In this approach, the waiting time for the next 
reaction to occur, as well as which specific reaction will occur 
are random variables that are completely determined by the mas- 
ter equation, so this approach should be considered the most ac- 
curate. In principle, multiple runs are needed to average out the 
random fluctuations, but in practice this is unnecessary if one 
only cares about the abundant species. This approach has been 
applied successfully t o astrochemical problems (Charnley 199*81 
1200 U IVasvunin et alJl2009l) even in the case of very large net- 
works |^^^uni^^l]l2002)- Besides providing results that are 
accurate, this approach is very easy to implement. However, it 
requires a very long run time for large networks if a long evolu- 
tion track is to be followed , although some approximate acceler- 
ated methods do exist (e.g. lGillespidl200^ . 

The SSA described above is somewhat different from a 
Monte Carlo (MC) approach which has al so been applied to 
astrochemistry (e.g. iTielens & Hagenlll982b : however, this ap- 
proach is not rigo rously consistent with the master equation (see 
the comment by ICharnlevI 2005h. and can l ead to a reaction 
probability higher than 1 (lAwad et al.l l2005l) in certain cases. 
The nomenclature of these two approaches is not always con- 
sistent i n the astrochem i cal ht erature^. For example, the SSA 
used by IVasvunin et alj (|2009|) is called the Monte Carlo ap- 
proach in their paper. Hereafter, we use the term "Monte Carlo" 
when referring to the rigorous stochastic simulation approach of 
Gillespie. 

By taking various moments of the maste r equation, the so- 
called moment equation (ME) is obtained dLips htat & BihanJ 
[2003; Barzel & Biham 2007a b). This set of equations describes 



the evolution of both the average population of each species and 
the average value of the products of the population of a group of 
species, usually cut off at the second order moments. Its formu- 
lation is similar to that of the RE, so it is relatively easy to im- 
plement. Furthermore, in this approach the gas phase chemistry 
and grain surface chemistry can be coupled together naturally. It 
has been tested on small surface networks. 

In the present paper, we propose yet another approach to 
modeling gas-grain chemistry, named the hybrid moment equa- 
tion (HME) approach. The goal is to find a systematic, auto- 
matic, and fast way to modeling gas-grain chemistry as accu- 
rately as possible. Our method is based on the ME approach. 
Different approximations are applied to the MEs at different time 
depending on the overall populations at that specific time. It 
is hybrid in the sense that the RE and the ME are combined 
together The_basic modification and competition scheme pre- 
sented in Garrodl(l2008l) can be viewed as a semi-steady-state ap- 
proximation to our approach (by assuming that the time deriva- 
tives of certain second order moments are equal to zero), while 
our approach can also be viewed as a combination of the ME ap- 
proach of Barzel & Biham (2007a) and the RE. In our approach, 
the MEs are generated automatically with the generating func- 
tion technique, and in principle MEs up to any order can be ob- 
tained this way. We benchmark our approach against the exact 
MC approach (i.e. the SSA of Gillespie). 

The remaining part of this paper is organized as follows. In 
section|2l we review the chemical master equation and ME, then 
describe the main steps of the HME approach. In section |3] we 
benchmark the HME approach with a cutoff at the second order 
and the RE approach against the MC approach with a large gas- 
grain network; we also tested the HME approach with a cutoff 
at the third order with a small network. In section we discuss 
the performance of the HME, and its relation with previous ap- 
proaches, as well as possibilities for additional improvements. 
Our way of generating the MEs is described in Appendix |A] 
A surface chemical network we used for benchmark is listed in 
AppendixlB] 



2. Description of the hybrid moment equation 
(HIVIE) approach 

In this section, we first review both the chemical master and mo- 
ment equations. Although t his content can be fou nd in many 
other papers (e.g., Chamlevi ll998t lGiIlesnidl2007h . we present 
them here as they are the basis of our HME approach. We then 
describe the MEs and REs for a simple set of reactions as an ex- 
ample, to demonstrate how the HME approach naturally arise as 
a combination of ME and RE. Finally we show the main steps of 
the HME approach. 

2.1. The chemical master equation and the moment 
equation (ME) 

A chemical system at a given time t can be described by a state 
vector X which changes with time, with its jth component Xj 
being the number of the jth species in this system. As a chemical 
system is usually stochastic, x should be viewed as a random 
variable, whose probability distribution function P{x, t) evolves 
with time according to the master equation (.Gillespie..2007i) 



For a discussion about the relations and di fferences between 
"stochastic simulation" and "Monte Carlo", see iKalos & WhTtiocg 
(120081) and lRiplevI i l2008h . 



d,P{x, t) = Y}ai{x - Vi)P{x - V,-, f) - ai(x)P(x)l 



(1) 



i=l 
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where a,(jr) is called the propensity function, a,(x)Af is the prob- 
ability that given a current state vector x an ith reaction will 
happen in the next infinitesimal time interval Af, and v, is the 
stoichiometry vector of the ith reaction. The sum is over all the 
reactions, and M is the total number of reactions. 

The ME is derived by taking moments of the master equa- 
tion. For example, for the^rsf order moment (xj), which is sim- 
ply the average number of species /', (xj ) = 2^. P(x, t)xj, its evo- 
lution is determined by (lGillespiell2007h 

d,{xj)^J]dAP(x,t)]xj 

X 

M 

= ^>["'(-*^ ~ Vi)P(x - Vi, t) - ai{x)P(x)] 

/ X 
M 

" Z Z [^""^ + t) - Xjai{x)P{x)] (2) 

/ X 

M M 

= ^ ^ vijai{x)P{x, = ^ Vijiuiix)), 

/ X ! 

where is the jth component of the stoichiometry vector of 
the ith reaction, i.e. the number of jth species produced (nega- 
tive when being consumed) by the ith reaction. For higher order 
moments, their coiTesponding evolution equations can be simi- 
larly derived, although the final form will be more complex. In 
AppendixlAl we present another method based on the generating 
function technique to derive the MEs, which is more suitable for 
programming. 

For the simplest network, in which all the reactions are 
single-body reactions, ai(x) is a linear function of x. In this case 
the ME is closed and can easily be solved. However, when two- 
body reactions are present, this is no longer true, as (fl,(jc)) might 
be of a form - 1)) or {x^xi), which is of order two and 

cannot be determined in general by the lower order moments. 
Hence additional equations governing their evolution should be 
included, i.e., they should be taken to be independent variables. 
The evolution equation of these second order moments may also 
involve moments of order three, and this process continues with- 
out an end, thus the ME is actually an infinite set of coupled 
equations (although in principle they are not completely inde- 
pendent if the chemical system being considered is finite, which 
leads to a finite-dimensional space of state vectors). The equa- 
tion cannot be solved without a compromise, e.g., a cutoff pro- 
cedure, except for the simplest cases in which an analytical so- 
lution is obtainable in the steady state. 

2.2. The MEs and REs for a set of reactions 

We take the following symbolic reactions as an illustrative ex- 
ample 



Adsorption: a 



hi 



Evaporation: A > 

Surface reaction: A + B 
Surface reaction: A+ A 



a, 

kAB 



C +D, 

E, 



(3) 
(4) 
(5) 
(6) 



where the ks are the reaction rates of each reaction, A - E are 
assumed to be surface species that are distinct from each other, 
and "a" is the gas phase counterpart of A. 



In the following we first write down the MEs and REs for this 
system, then discuss the relations and differences between them, 
as well as the relation between a cutoff of MEs and a cutoff of 
master equations in previous studies. These discussions will be 
essential to developing our HME approach. 

2.2.1 . The MEs for this system 

The propensity functions for the above four reactions are k^dci, 
^evap'^, kABAB, and k^AMA - 1), respectively. Here for conve- 
nience we use the letter "A" to represent both the name of a 
species and the population of the corresponding species. 
For the first order moments, we have 

d,{A) =^,d<fl) - ^evap<A> - kAB{AB) - 2A:aa<A(A - 1)>, (7) 
d,{C) =^ab(AB>, (8) 
d,{E) =^AA<A(A - D). (9) 

Other similar equations are omitted. The symbol (*) is used to 
represent the average population of "*" in the system; the aver- 
age should be understood as an ensemble average. The second 
order moments (AB) and (A (A - 1)) have their own evolution 
equations, which are 



d,{AB} =k.,i{aB} - ^evap<AB) 

-A:ab[<A(A-1)B) + <AB(B-1))- 
- 2^AA<A(A - 1)B), 



(AB)] (10) 



dr{A{A - 1)) =2^,d(flA) - 2feevap<A(A - 1)) 

- 2^ab(A(A - 1)B) (11) 

- 2^aa[2<A(A - 1)(A - 2)) + (A(A - 1))]. 

For this simple example set of reactions (equation (l3]-|6) ), the 
above equatio ns can be easily obtain ed from the master equa- 
tion (see, e.g.. iLipshtat & BiharDl l2003. page 8). In the general 
case (e.g., when A - E are not completely distinct from each 
other), an automatic way of obtaining the MEs is described in 
Appendix [a] The method described there is also applicable to 
moments with any order, and to all the common reaction types 
in astrochemistry. 

In general, the third order moments in the above equations 
cannot be expressed as a function of the lower order moments, so 
they need their own differential equations. In the case of a cutoff 
at the second order, the chain of equations, however, stops here. 
We describe the method required to evaluate them in section 1231 



2.2.2. The REs for this system 

When using REs, equations (l7l-[TTI) are replaced by 

dr{A) = k.,d{a) - ^evap<A) - ^AB<A)(B) - 2^AA<A>^ 0) 

d,{C) = A:ab<A)<B>, (EI) 

dr{E) = kAAiAf. m 

d,[{A){B)] = k.,d{a){B) - ^evap<A)(B) 

- kAB[{Af{B) + {A){Bf] m) 



- 2kAA{Af{B), 
d,[{Af] = 2k,A{a){A) - 2A:evap<A>2 

-2^ab<A>2(B>-4^aa<A>^ 



(HU) 



The equations for {A){B) and (A)^ are of course not needed in the 
RE approach but are simply derived from equation (|7]) (and an 
omitted similar equation for {B)) using the chain rule of calculus. 
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2.2.3. The relation between MEs and REs 

The differences between the MEs (equationl7]-[TTb and the REs 
(equation |7] -[H]) in the present case are as follows: All the 
(AB) are replaced by {A){B), all the {A(A - 1)) are replaced by 
<A)2, all the {A(A - l)B) are replaced by {Af{B), the {AB(B - 1)) 
is replaced by {A){Bf-, and the (A(A - 1)(A - 2)) is replaced by 
{A)^ . Furthermore, the term kp^^iAB) in equation ( fTOt and the 
term k/^/^{A{A - 1)) in equation (fTIT i disappear in the RE (fTOi ) 
and dni). 

These differences make clear why the REs are accurate when 
the involved species are abundant (namely when (A)»l and 
This is because, in this case, {AB) can be approximated 

well 

b>EI {A){B), and (A(A - 1)) can be approximated well by 

{A)\ 

The RE approach will be erroneous when (A) or {B) are 
smaller than 1 because, in this case, the correlation between A 
and B might cause {AB) to differ considerably from {A){B), and 
the fluctuation in A might cause (A (A - 1)) to differ consider- 
ably from (A)^. It can also be viewed like this: in equation ( fTOi ) 
and equation ( fTTI ) that govern the evolution of second order mo- 
ments, the omitted term A:ab(AB> might be much larger than the 
retained terms such as (A)^(B) or (A)(B)^, and the omitted term 
kp,/^{A{A - 1)> might be much larger than the retained term (A)^. 

2.2.4. The relation between a cutoff of MEs and a cutoff of 
possible states in previous master equation 
approaches 

In Eqs. dTl-fTTI) we do not write terms such as (A(A - 1)> in the 
split form (A^) - (A). We keep terms such as (A(A - \)B) and 
(A(A - 1)(A -2)) in their present forms intentionally. One reason 
for this is that terms such as (A (A - 1)) look more succinct and 
follow naturally from our way of deriving them (see Appendix 
lAt . When (A) » 1, (A(A - 1)) and {AB) can be directly replaced 
by (A)^ and {A){B), respectively, to obtain the RE formulation. 

More importantly, this formulation can be directly connected 
to the cutoff schemes in the previous master equa tion approaches 
(e.g., lBihametalJl200II: IStantcheva e"tal1l2002l) . For example, 
in a scheme in which no more than two particles of A are ex- 
pected to be present on a single grain at the same time, we have 
P(A>2) = 0. In this case, <A(A - 1)(A - 2)> = 21=3 P(A)A(A - 
1)(A - 2) = 0. Thus we see that a cutoff at a population of two in 
the master equation approach corresponds naturally to assigning 
a zero value to moments containing A more than twice, as far as 
the moments are defined in the form presented above. 

2.3. The HME approach 

The HME approach is a combination of the ME and RE ap- 
proaches. The basic idea is that, for deterministic (average pop- 
ulation >1) species, the REs are used, while for stochastic (av- 
erage population <1) species, the stochastic effects are taken 
into account by including higher order moments in the equa- 
tions. Since a deterministic species may become stochastic as 
time goes by, and vice versa, the set of ODEs governing the evo- 
lution of the system also changes with time, and is determined 
dynamically. A flow chart of our HME code is shown in Fig.[T] 

^ Assuming Poisson statistics, we have 

\{AB) - {A){B)\ ^ l_L + _L « 1 
{A){B) ~ V <A> <S> 



We first set up all potentially needed MEs (using the proce- 
dure described in Appendix |A|l, with a cutoff of moments at a 
prescribed order (usually two). After this and some other initial- 
ization work, the program enters the main loop. 

The main loop contains an ODE solver because the system 
of MEs is a set of ordinary differential equations (ODEs). We 
use the solver from the ODEPACK packag^ 

Not all MEs and moments are used at all times; which ones 
are used is determined dynamically. In each iteration of the main 
loop, we verify whether some surface species have changed from 
stochastic to deterministic, or from deterministic to stochastic. 
The gas phase species are always treated as deterministic, re- 
gardless of how small their average populations are. In either of 
these two cases, we re-examine all the moments, and determine 
the way to treat them. There are four cases: 

1. All the first order moments are treated as independent vari- 
ables. 

2. If a moment consists of only stochastic species, and its order 
is no larger than the prescribed highest allowed order, it will 
be treated as an independent variable, and the corresponding 
moment equation will be included and solved. For the sake 
of numerical stability, its value should be no larger than its 
deterministic counterpart. For example, if the ODE solver 
yields a value of {AB) > {A){B), then the latter value will be 
assigned to {AB). 

3. If a moment consists of only stochastic species, and its order 
is larger than the prescribed highest allowed order, its value 
will be set to zero, and of course, its moment equation will 
not be solved. For example, if (A)<1 and {B)<1, then, with 
a highest allowed order set to two, moments such as (A(A - 
1)(A - 2)) and <A(A - l)B) w ill be set to zero. This follows 
from the discussion in section l2".2.4l 

4. If a moment contains at least one deterministic species, it 
will not be treated as an independent variable, and its mo- 
ment equation will not be solved. It can be evaluated in 
the following way: assuming that the moment under con- 
sideration has a form {AB(B - 1)), and that A is determin- 
istic (i.e. (A)>1), then the value of {AB{B - 1)> is set to be 
{A){B(B - 1)). If B is also deterministic, then it will be eval- 
uated as {A){B)^. This follows from the discussion in section 
IZ23] 

From these procedures, we see that the number of equations, as 
well as the form of these equations will change when a transi- 
tion between stochastic and deterministic state of certain species 
occurs. Each time the ODE system is updated, the ODE solver 
must therefore be re-initialized. 

It seems possible to replace the sharp transition between the 
stochastic and deterministic state of a species (based on whether 
its average population is smaller than 1 ) with a s mooth transition , 
e.g., using a weight function similar to that in iGarrodI (l2008h . 
However, it is not mathematically clear which weight function 
we should choose, and an arbitrary one might cause some artifi- 
cial effects, so we prefer not to use this formulation. 

3. Benchmark with the Monte Carlo approach 

We compare the results of our HME approach with those from 
the exact stochastic simulation ( Gillespidl2007l ICharnlevll 19981 
I2OOI; .Vasyunin et al. 2009). The RE results are also com- 
pared for reference. As in previous studies (IChamlevI 1200 1[ 
IVasvunin et al.ll2009l) . we consider a closed chemical system in 

* Downloaded from www.netlib.org 
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1. Import all the physical parameters 
and reactions, calculate the rates. 



2. Prepare all the potentially needed 
moment equations (see Appendix A). 



3. Initialize all the moments, assuming no 
species is on the grain at the beginning. 



C 



4. Start iteration. 



J 




5. Stochasticity changed? 
|yes 



6. Determine which moments to 
be kept and which to be dropped. 



( 7. Re-initialize the ODE solver. ] 



c 



T 



8. CaU the ODE solver. 



T 



> 




9. Final time reached? 
I yes 



10. Finish. 



Fig. 1. A flow chart of the main components of our HME code. 
Steps 2, 5, 6, and 8 are described in detail in the text. 



a volume containing exactly one grain particle. The number of 
each species in this volume is called a "population", which can 
be translated into an abundance relative to H nuclei by multiply- 
ing it by the dust-to-gas ratio, which is 2.8 x 10~'^(0.1 /im/r)-', 
where r is the grain radius, assuming an average molecular 
weight of 1.4, a dust-to-gas mass ratio 0.01, and a density of 
grain material of 2 g cm"-'. 

In the MC approach, the number of each species in this vol- 
ume at any time is an integer Owing to the large amount of 
time steps (>10^), it is impractical to store all intermediate steps, 
so we average the population of each species in time, weighted 
by the time intervals (remember that the lengths of time inter- 
vals between reactions are also random in MC). Because of this 
weighted average (rather than merely saving the state vector at 
certain instants), the MC approach can resolve average popula- 
tions much smaller than one, although the fluctuations that are 
intrinsic to the MC approach can be larger than the average pop- 
ulations when the latter is small. 

We first demonstrate how our method works for a large gas- 
grain network. We then show that for a small surface network, 
third order moments can also be included to improve the accu- 
racy. 



3.1. Test of the HME approach truncated at the second order 
on a large gas-grain network 

We use the "dipole-enhanced" version of the RATE06 gas phase 
reaction networ^ ( Woodall et al. 2007), coupled with a surface 
network of iKeanel d 1 9971) (see Appendix iBl). The surface net- 
work contains 44 reactions between 43 species, which is ba- 
sic ally a reduced and slig htly revised version of the network 
of ITielens & HagenI (Il982h . containing the formation routes of 
the most common grain mantle species, such as H2O, CH3OH, 
CH4, NH3, etc. This surface network is not really large in com- 
pa rison with some of the previous works, such as that used 
by iGarrod et al.l (l2009h . However, it is already essential for the 
most important species. The energy barriers for thermal desorp- 
tion and diffusion are taken from Stantcheva & Herb st {2003- 
Diff'usion of H atoms on the surface through quantum tunneling 
is included. Desorption by cosmic rays is taken into account fol- 
lowing the approach of Hasegawa & Herbst (1993). The rate co- 
effic ients of the gas p hase reactions are calculated according to 
Woodall et al. while the rate coefficients of the surface 

(119921). 



Hasegawa et al 



The 



reactions are calculated followin _ . 

initial condition is the same as in lStantcheva & HerbstI (l2004l) 

We assumed a dust-to-gas mass ratio of 0.0 1 . The grain mass 
density is taken to be 2 g cm"^, with a site density 5 x lO'^ cm~^. 
Two grain sizes have been used: 0. 1 fim and 0.02 yum. A cosmic 
ray ionization rate of 1.3 x 10"'^ s"' is adopted. Four different 
temperatures (10, 20, 30, 50 K) and three different densities (10^, 
10"*, 10^ cm"^) have been used. In total, the comparison has been 
made for 24 different sets of physical parameters. These con- 
ditions are commonly seen in translucent clouds and cold dark 

clouds. 

As in lGarrod et al.l (12009). we make a global comparison be- 
tween the results of MC, HME, and RE. For each set of physical 
parameters, the comparisons are made at a time of 10^, lO'*, and 
10^ years. We calculate the percentage of species for which the 
agreement between MC and HME/RE is within a factor of 2 
or 10. Only species with a population (either from MC or from 
HME/RE) larger than 10 are included for comparison. This is 
because for species with smaller populations, the intrinsic fluc- 
tuation in the MC results can be significant. For several different 
sets of physical parameters, we repeated the MC several times to 
get a feeling for how large the fluctuation magnitude would be, 
although this is impractical for all the cases. 

The comparison results are shown in Table [1] (grain radius 
= 0.1 fim) and Table |2] (grain radius = 0.02 jjm). The HME ap- 
proach always has a better global agreement (or the same for 
several cases) than the RE approach in the cases we tested. The 
typical time evolution of certain species is shown in Fig. |2] In 
each panel of the figure, the species with a name preceded with 
a "g" means it is a surface species. 

The poorest agreement of HME (Fig. O is at f = 10^ year 
for T=20 K, riu = 10^ cm~^, and grain radius - 0.02 fim. This 
is mostly because at the time of comparison the populations of 
certain species were changing very rapidly, so a slight mismatch 
in time leads to a large discrepancy. This mismatch is probably 
caused by the truncation of higher-order terms in the HME (see 
section [XTb . For gN2 in Fig. |3] its population seems to be sys- 
tematically smaller in HME than in MC during the early period, 
although the HME result matches the one from MC at a later 
stage (after 3 x lO-' years). 

The RE is as effective as the HME in several cases, when the 
temperature is either relatively low (~10 K) or high (~ 50 K) 



' http://www.udfa.net/ 
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(see also IVasvunin et af]l2009l) . and generally works better for a 
grain radius of 0. 1 yum than of 0.02 jjm. When the temperature is 
very low, many surface reactions with barriers cannot happen (at 
least in the considered timescales). On the other hand, when the 
temperature is high, the surface species evaporate very quickly 
and the surface reactions are also unable to occur. In these two 
extreme cases, the surface processes are inactive, and the RE 
works fine. 

The RE becomes problematic in the intermediate cases, 
when the temperature is high enough for many surface reactions 
to occur, but not too high to evaporate all the surface species; in 
these cases the HME represents a major improvement over the 
RE. For a smaller grain radius, the population of each species in 
a volume containing one grain will be smaller, thus the stochas- 
tic effect will play a more important role, and the RE will tend 
to fail. 

We note that, in the HME approach, there is no elemental 
leakage except those caused by the finite precision of the com- 
puter In all the models that we have run, all the elements (includ- 
ing electric charge) are conserved with a relative error smaller 
than SxlO"'**. The reason why elemental conservation is always 
guaranteed is that either the rate equations or the moment equa- 
tions for the first order moments conserve the elements. 

When comparing the results from the HME approach with 
those from MC simulation, it is important to see how the in- 
trinsic fluctuation in MC behaves. If we assume the probability 
distribution of the population of a species, say A, is Poissonian, 
then the variance of A is cr^(A) = (A). Hence if (A) is small, 
the relative fluctuation of the MC result can be quite large. This 
fluctuation might be smoothed out by means of a weighted av- 
erage in time, but this procedure does not always work. This is 
why we choose to only compare species with a population higher 
than 10, corresponding to an abundance relative to H nuclei of 
2.8x10"" (for grain radius =0.1 /im) or 3.5 x 10"^ (for grain ra- 
dius = 0.02 iim). For a real reaction network, it is usually diflicult 
to predict the intrinsic fluctuation in a MC simulation, unless it 
is repeated many times. These fluctuations will not have any ob- 
servational eff'ects, because along a line of sight there are always 
a large number of a certain species (as far as it is detectable) and 
the fluctuations are averaged out. 

We note that the gas phase processes are not treated iden- 
tically in our HME approach and MC simulation. In the MC 
approach, the gas phase processes a re always treated as being 
stochastic (see, e.g., Charnlev 1998; V asvunin e t al. 2009), in the 
same way as the surface processes. However, in our HME ap- 
proach, the gas phase species are treated in a deterministic way, 
i.e., REs are always applied to them. This means that even if two 
reacting gas phase species A and B both have average popula- 
tions much smaller than one, we still assume that (AB) - {A){B). 
This is physically quite reasonable, because the presence of large 
amounts of reacting partners in the gas phase (if n ot limited 
to a vo lume containing only one dust grain; see, e.g.. lChamlevl 
ensures that the RE is applicable. However, although it 
might sound a bit pedantic, mathematically this is not equivalent 
to the MC approach, and some discrepancies caused by this are 
expected. For a large network, it is impractical to treat the gas 
phase processes in the same way as the surface processes in the 
HME approach, because in that case the number of independent 
variables in the ODE system will be quite large (at least no less 
than the number of two-body reactions), and the performance of 
the ODE solver will be degraded. 



3.2. Test of the HME approach truncated at the third order 
on a small surface network 

To test the improvement in accuracy when the cutoff is made at a 
higher order, we compare the results of the HME approach with 
a cutoff at the second order to those obtained from the same ap- 
proach with a cutof f at the third order We use a small surface 
reaction network of IStantcheva & Herbsd (l2004l) . containing 17 
surface reactions between 21 species, producing H2O, CH3OH, 
CH4, NH3, and CO2. No gas phase reactions are included, ex- 
cept adsorption and desorption processes. The initial gas phase 
abundances of the relevant species are obtained from the steady 
state solution of the RATE06 network under the corresponding 
physical conditions. 

As before, we run the HME, RE, and the MC code for diff'er- 
ent sets of physical parameters. Although by transferring from 
the RE to the second order HME a major improvement in accu- 
racy can be obtained, the inclusion of the third order moments to 
the HME usually only improves slightly over the second order 
case. In Fig. HI we show an example (T = 10 K, hh = 2 x 10^ 
cm"^, grain radius=0.02 fim), in which the distinctions between 
the results from the second and third order HME are relatively 
large. 

For several species, we note that the third order HME is still 
unable to match the MC results perfectly, and for gHCO (Fig.|4|i 
the third order HME even produces an artificial spike in the time 
evolution curve. The results from the third order HME are oth- 
erwise of greater accuracy than the second order one, the abun- 
dances of gHiCO and gCH30H in particular being in almost 
perfect agreement with those from the MC approach. In the case 
of gHCO, the timescale mismatch between HME and MC is al- 
leviated by including the third order moments. 

It might be useful to see the difference between the second 
order HME and the third order HME in a computational sense. 
For the current reaction network with physical parameters de- 
scribed above, the number of variables (same as the number of 
equations, which changes with time) is 145 initially in the sec- 
ond order case, and this number becomes 705 for the third order 
case. To reach a time span of ~10^ year, the second order HME 
takes about 3 seconds, while the third order one takes about 220 
seconds on a standard desktop computer (a CPU @ 3.00 GHz 
with double cores, 4 GBytes memory). The number of variables 
depends on the network structure, and it is not straightforward 
to derive a formula to calculate it. Qualitatively, this number (as 
a function of the number of reactions or the number of species) 
seems to increase with the cutoff order less quickly than expo- 
nential growth. However, such a "mild" increase affects the be- 
havior of the ODE solver quite significantly. This is partly be- 
cause the solver contains operations (such as matrix inversion) 
that become slower as the number of variables become larger 
An increase in the number of variables might also increase the 
stiffness of the problem, let alone the memory limitation of the 
computer. For the larger reaction network described in the pre- 
vious section, the third order HME would involve about 5000 
variables and has not been tested successfully. 

4. Discussion 

In our HME approach, we have used a general and automatic 
way to derive the MEs. For large gas-grain networks, the MEs 
are cut off at the second order For small networks, a cutoff at the 
third order is possible and higher accuracy can be obtained. We 
incorporate a switching scheme between the ME and RE when 
the average population of a species reaches 1 . 
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Table 1. Percentage of agreement between the results from MC and those from HME/RE. The comparison is only made between 
those species with populations (from MC or HME/RE) larger than 10. The two numbers in each table entry means the percentage 
of agreement within a factor of 2 or 10, respectively. The grain radius is taken to be 0.1 fim. 
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Table 2. Same as Table [T| except a smaller grain radius of 0.02 fim is taken. 
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Fig. 2. Typical time evolution of the average populations of certain species from MC (solid lines), HME to the 2 order (dotted 
lines), and RE (dashed lines). Note that the Monte Carlo has been repeated twice. The y-axis is the number of each species in a 
volume containing exactly one grain. To translate it into abundance relative to H nuclei, it should be multiplied by 2.8 x 10 
Physical parameters used: T = 20 K, n = 2 x 10^ cm"-', grain radius = 0.1 yum. 
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10' 10^ 10* 10" 
Time (year) 



•Q^ to* 
Time (year) 



Fig. 3. Cases in which the agreement between the results of MC and those of HME are not so good, especially at t = lO' year. The 
y-axis is the number of each species in a volume containing exactly one grain. To translate it into abundance relative to H nuclei, it 
should be multiplied by 3.5 x 10 Physical parameters used: T = 20 K, n = 2 x 10^ cm"^^, grain radius = 0.02 fj.m. 



The results from HME are more accurate than those from 
the RE in the cases we have tested, when benchmarked against 
the exact MC results. The abundances of almost all the abundant 
species (>2.8xl0"" for a grain radius of 0.1 jum and >3.5xlO"^ 
for a grain radius of 0.02 yum) from HME are accurate to within 
a factor of two, especially at later stages of the chemical evolu- 
tion, while in some cases nearly 40% of the results from RE are 
incorrect by a factor of at least ten. 

In terms of computation time, our approach usually takes 
several tens of minutes to reach a evolution time span of 10*' 
years, so it is slower than the RE, but faster than the MC ap- 
proach (which usually takes from several hours to day s). Our 
approa ch may also be slower than the MRE approach of lGarrodI 
(I2008h . because more variables (namely the moments with or- 
ders higher than one) are present in our method, and the ODE 
system in HME is usually stiffer. For example, for a moderate 
temperature many surface reactions can be much faster than any 
gas phase reactions, and yield a very large coefficient in some of 
the MEs. However, this is not the case in the stochastic regime of 
the MRE approach, because when a competition scheme is used 
in MRE, such a large coefficient does not appear. In t his sense, 
we also advocate the MRE approach of iGarroJ (l2008h . 

Mathematically, our app roach is partially e quivalent to the 
master equation approach of IStantcheva & Herb st (2004) in two 
respects. 1) They separated stochastic and deterministic species, 
which is similar to our adopting RE for the abundant species. 2) 
They set a cutoff for the possible states of the stochastic species. 
This is in essence equivalent to letting moments containing these 
species with order higher than a certain number equal to zero . 

Our appro ach can also be vie wed as a combination of some 
of the ideas of lGarrodI (l2008l) and lBarzel & BihamI (l2007al) . The 



basic modification and competition scheme in lGarrodI (1200 8') can 
be derived from the MEs, wit h a semi-steady s tate as sumption 
for the second order moments. iBarzel & BihamI (l2007ah used the 
MEs, but they did not include a switch scheme, and their way of 
deriving the MEs is different from the one in the present paper. 

There are still many possibilities for improvement. Although 
in principle moments with any order can be included, the num- 
ber of equations grows quite quickly with the cutoff order, which 
makes the system of equations intractable with a normal desk- 
top computer. It is unclear whether it is possible to include the 
moments selectively. It is unclear whether there are better, and 
more mathematically well founded strategies than the switch at 
an average population of 1 . The present approach is usually sta- 
ble numerically. However, this is not always guaranteed, espe- 
cially if higher order moments are to be included. The behavior 
of the numerical solution also depends on other factors, such as 
the ODE solver being used and the tunable parameters for it, 
while the MC approach does not have such issues. In this sense, 
the MC approach is the most robust. 

Even in the accurate MC approach described above, the de- 
tailed morphology of the grain surface and the detailed reac- 
tion mechanism is not taken into account. One step in this di- 
rection would be to take into accou nt the lay e red st ructure of 
the grain mantle. Th i s was done by ICharnlevI (l200lh (see also 
IChamlev & RodgersI (|2009|) ) by means of stochastic simulation. 
It could also be included in the HME approach, as far as the 
underlying physical mechanism could be described by a master 
equation. 

However, a microscopic MC approach has also been used 
to study the grain ch emistry (see, e.g., IChang et al] 120051: 
ICuppen & Herbstl2007l) . In this approach, the morphology of the 
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Fig. 4. Comparison of the results from MC, HME to the 2°'' order, HME to the 3"^ order, and RE. The y-axis is the number of each 
species in a volume containing exactly one grain. To translate it into abundance relative to H nuclei, it should be multiplied by 
3.5 X 10 Physical parameters used when running these models include T = 10 K, «h = 2 x 10^ cm grain radius - 0.02 fim. 



grain mantle and the interaction between species are modeled in 
detail. As far as we know, this approach is only practical when 
the network is small. It remains unclear whether it is possible to 
incorporate these details into the current HME approach. 

In some cases, errors caused by uncertainties in the reaction 
mechanism and rate parameters might be larger than those intro- 
duced by the modeling method (Vasvuninet al. 2008). Hence, 
further experimental study and a more sophisticated way of in- 
terpreting those results would be indispensable. 
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Appendix A: A method to generate the moment 
equations based on the generating function 

We describe our means of getting the MEs. Our method is au- 
tomatic, and can be easily coded into a computer program. It 
is applicable to moments of any order and all the common as- 
trochemical reactions. It makes use of the probability generat- 
ing function. While pr eparing the present paper, we noted that 
iBarzel & BihamI 42010) also proposed a binomial formulation 
of ME, which in essence is partly equivalent to our approach 
presented here, although our way of deriving the MEs is quite 
different from theirs. 

For a probability distribut ion P(x, t), the corre sponding gen- 
erating function is defined as dvan Kampenll2007h 



f(ZuZ2,...,t) = J]Pix,t)zfzf-- 



(A.l) 



Here all the z,s should be thought of as merely symbols without 
any physical meaning, and they have a one-to-one correspon- 
dence with the x,s. 

It is obvious that f(z - 1, f) = 1, which is just the normal- 
ization condition for probability. It is also easy to see that the 
average population of the ith species is 



(Xi) = d:.f(zl,Z2,---,t)\z=l- 



(A.2) 



The right hand side of the above equation means taking the par- 
tial derivative first, then assigning a value one to all the z,s. 

For the second order moment between two distinct species / 
and j, we have 

{xiXj)^d,,d-J(zuZ2,...,t)\z=i- (A.3) 
If i equals / in the above equation, then what we actually get is 

{xi(xi-l))^dlf(zuZ2,.. .,1)1=1. (A.4) 



In general, we have 

{XiXjXk ■ ■ ■ > = 



■/(zi,Z2,...,f)lz=i- 



(A.5) 



If several of the subscripts are the same in the left hand side of 
the above equation, say, i - j - k, then the second should be 
understood as (xj - 1), while the third should be understood as 
(xii - 2), and so on. 

From the master equation in Eq. ([TJ, it seems possible to 
get an equation for the evolution of f(z, t) in the general case. 
However, if the propensity functions aiix) (see Eq. ([TJ) are al- 
lowed to take any functional form, then this is not straightfor- 
ward. Fortunately, in practice ai{x) usually has a very simple 
form. On the other hand, we note that in the right hand side of 
the master equation in Eq. ([T]i the contributions from all the reac- 
tions are added linearly. Hence the contribution of each reaction 
to the evolution of generating function can be considered inde- 
pendently of each other. 

We assume there is only one reaction in the network, which 
has a form 



should be understood as explained in the sentence following Eq. 
(IA.5I )). then the generating functiorQ will evolve according to 



dtf = k{yiy2 ■■■ym- xiX2 ■ ■ ■ Xn)d^,d_.,, ■ ■ ■ d^^J. 



{A.l) 



It is not diflicult to derive the above equation from the master 
equation and the definition of generating equation and our as- 
sumption about the propensity function. 

We note that equation ( IA.7l i has a very simple pattern that is 
easy to remember: (a) The constant coefficient is the rate coeffi- 
cient; (b) In the parenthesis, the symbols of all the products are 
multiplied together, with a coefficient -i-l, while the symbols of 
all the reactants are multiplied together, with a coefficient -1; (c) 
In the differential part, all the reactants are present as they are in 
the left hand side of equation (IA.6I 1, while none of the products 
appear. 

We now attempt to derive the ME. We obtain the evolution 
equation of each moment (as defined in Eq. dA.Sb ) by simply 
differentiating both sides of equation (IA.7b with respect to the 
relevant components in the moment, then setting all the symbols 
to a value of one. 



For example, for the reaction A + B ■ 
equation of the generating function / is 



C + D, the evolution 



d,f = k(CD - AB)dAdBf. 

For {AB), we differentiate both sides of the above equation by A 
and B. We obtain 

d,dAdBf ^k[{CD-AB)d\dlf 

- Adidef - BdAdlf - dAdefl 

Next we assign a value of one to all the symbols (A - D) appear- 
ing in the resulting expressions, and "translate" the remaining 
terms into moments (recalling the remark about equation (IA.5b ). 
obtaining 

d,{AB) = -k[{A(A - 1)B) + {ABiB - 1)) + {AB)]. 

Although the above derivation involves differentiations, 
these operations can be easily translated into some combinatorial 
rules and written as a computer program. A recursive procedure 
is needed to generate all the potentially needed moments up to a 
given order. 



X\ + X2 + ■ 



■yi+y2 + ---+yn. 



(A.6) 



where the x,s and y,s represent the reactants and products, 
which do not have to be different from each other. We also use 
these symbols to represent the populations of the coiTesponding 
species. If, given a population of xi, jca, . . . , x„, the probability 
that the above reaction will happen in a unit time is kx\X2 ■ ■ ■ x„ 
(namely, the propensity function a{x) - kx\X2 ■ ■ ■ x„; the product 



* Instead of using z, s as symbols for the independent variables of the 
generating function /, we use x,s and y,s instead. This will not cause 
any confusion. 
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Appendix B: A surface reaction networi< we used to 
test our code 



Table B.l. The surface network used in Section lTTI of this paper. 
The original references listed below this table should be cited if 
these data are to be used. We note that the validity of the numer- 
ical method (namely the HME approach) presented in this paper 
does not depend on the specific test network that we used. 



Number 


Reactants 


Products 




£.eac (K) 


1 


H 


H 


H2 




0.0 


2 


H 


o 


OH 




0.0 


3 


H 


OH 


H2O 




0.0 


4 


H 


On 


O2H 




1200.0 


5 


H 


OoH 


H2O2 




0.0 


6 


H 


Ht02 


H2O 


OH 


1400.0 


7 


H 


Oj 


O2 


OH 


450.0 


8 


H 


CO 


HCO 




1000.0 


9 


H 


HCO 


H2CO 




0.0 


10 


H 


HiCO 


CH3O 




1500.0 


11 


H 


H2CO 


H2COH 




1500.0 


12 


H 


CH^iO 


CH,OH 




0.0 


13 


H 


HtCOH 


CH,OH 




0.0 


14 


H 


HCOO 


HCOOH 




0.0 


15 


H 


N 


NH 




0.0 


16 


H 


NH 


NH, 




0.0 


17 


H 


NH7 


NH3 




0.0 


18 


H 


c 


CH 




0.0 


19 


H 


CH 


CHo 




0.0 


20 


H 


CHt 


CH, 




0.0 


21 


H 


CH3 


CH4 




0.0 


22 


H 


CN 


HCN 




0.0 


23 


H 


NO 


HNO 




0.0 


24 


H 


NOt 


HNO2 




0.0 


25 


H 


NO3 


HNO, 




0.0 


26 


H 


NHCO 


NH2CO 




0.0 


27 


H 


NH2CO 


NH2CHO 




0.0 


28 


H 


NjH 


N2H2 




0.0 


29 


H 


N2H2 


NjH 


H2 


650.0 


30 








O7 




0.0 


31 





O2 


O3 




1200.0 


32 





CO 


CO2 




1000.0 


33 


O 


HCO 


HCOO 




0.0 


34 





N 


NO 




0.0 


35 





NO 


NO2 




0.0 


36 





NO2 


NO3 




0.0 


37 


O 


CN 


OCN 




0.0 


38 


C 


N 


CN 




0.0 


39 


N 


N 


N2 




0.0 


40 


N 


NH 


NjH 




0.0 


41 


N 


HCO 


NHCO 




0.0 


42 


H2 


OH 


H2O 


H 


2600.0 


43 





HCO 


CO2 


H 


0.0 


44 


OH 


CO 


CO2 


H 


80.0 



References. iKeand < 19971) and lTielens & Haged <1982l') . 



